C.................... F107Kp-Update.for .......................
C..... This program reads the IRI APF10.7 index file and outputs the FLIP F10Kp file.
C..... This program was originally written to read the indices file from NGDC but
C..... NOAA discontinued the process on May 1, 2018.
C..... Note that some IRI daily F10.7 values are bad and replaced with NGDC values
C..... Also some 81-day values are bad. Program recalculates all values, except the 
C..... first 41 and the last 41.
C..... Rewritten in May 2019 by Phil Richards
      IMPLICIT NONE
      INTEGER I,IDIM,J,K                  !.. loop control variables
      PARAMETER (IDIM=55000)              !.. Array dimension. Big enough for 150 years
      INTEGER YEAR(IDIM),MONTH(IDIM),DAY(IDIM)
      INTEGER YYYY(IDIM)                  !.. Four digit year array         
      INTEGER YADD                        !.. Used to convert from 2 digit to 4 digit year
      INTEGER DDD                         !.. Day of the year to calculate Sun-Earth distance
      INTEGER Kp(IDIM,9),Ap(9)            !.. Kp, Ap values
      INTEGER APtoKP(0:400)               !.. Ap to Kp conversion values
      REAL F107YR                         !.. Yearly average F10.7
      INTEGER QB                          !.. Dummy for read
      REAL F107ADJ(IDIM),F107ADJA         !.. Adjusted F10.7 indices. No longer used
      REAL F107OBS(IDIM)                  !.. Observed F10.7 index observed
      REAL F107AVE,F107OBSB               !.. F10.7 81-day averages (IRI,PGR)
      REAL SEDIST(IDIM)                   !.. Sun-Earth distance factor
      REAL F107OBSA(IDIM)                 !.. 81 day average of observed F10.7

      !.. Ap to Kp conversion
      DATA APtoKp/2*0,3,7,10,13,17,2*20,3*23,3*27,3*30,
     >  4*33,5*37,5*40,7*43,9*47,8*50,11*53,13*57,14*60,17*63,
     >  21*67,22*70,25*73,28*77,29*80,64*83,100*87,90/
      
      CALL OPEN_FILE()      !... open files assigned in batch (.bat) file

      DDD=1      !.. Keeps track of the day of year
      YADD=1900  !.. Convert year from two to four digit
      K=1        !.. Counts the number of good data lines read in

      !.. Priming Read to skip any headers and get first line of the IRI data file
 10   READ(1,'(13I3,3F5.1)',END=91,ERR=10) YEAR(K), MONTH(K), DAY(K),
     >   (AP(I),I=1,9),QB,F107ADJ(K),F107AVE,F107YR

       
      !.. Loop to process the lines of data
      DO J=2,IDIM
         DO I=1,9
             KP(K,I)=APtoKP(AP(I))   !.. convert AP to KP
         ENDDO

         !... Covert to 4 digit year         
         IF(YEAR(K).LT.YEAR(K-1)) YADD=2000  !.. Change to 2000s when YY goes from 99 to 00
         IF(YEAR(K).NE.YEAR(K-1)) DDD=1      !.. Reset day number when year changes
         YYYY(K)=YEAR(K)+YADD            !.. Convert from 2 digit to 4 digit year

         !.. Take care of known bad data up to 2019 using NGDC values
         IF(YYYY(K).EQ.2005.AND.MONTH(K).EQ.8.AND.DAY(K).EQ.22) 
     >      F107ADJ(K)=107.22
         IF(YYYY(K).EQ.2005.AND.MONTH(K).EQ.9.AND.DAY(K).EQ.9) 
     >      F107ADJ(K)=100.62
         IF(YYYY(K).EQ.2005.AND.MONTH(K).EQ.9.AND.DAY(K).EQ.13) 
     >      F107ADJ(K)=115.02
         IF(YYYY(K).EQ.2006.AND.MONTH(K).EQ.12.AND.DAY(K).EQ.6) 
     >      F107ADJ(K)=99.72
         IF(YYYY(K).EQ.2011.AND.MONTH(K).EQ.3.AND.DAY(K).EQ.7) 
     >      F107ADJ(K)=150.40
         !.. Test to see that the year month and day are continuous
         IF(K.GT.1) CALL TEST_DATA(IDIM,K,YYYY,MONTH,DAY)

         !.. test F10.7 is within reason
         CALL TEST_F107(IDIM,K,YYYY,MONTH,DAY,F107ADJ)

         !.. factor to adjust F10.7 for the Sun-Earth distance
         SEDIST(K)=(1.0+0.0167*COS(6.283185*REAL(DDD-3)/365.25))/
     >      (1.0-0.0167*0.0167)

         F107OBS(K)=F107ADJ(K)*SEDIST(K)**2  !.. IRI observed F10.7
         F107OBSA(K)=F107AVE*SEDIST(K)**2    !.. IRI 81-day average observed F10.7

         K=K+1       !.. increase the number of lines processed
         DDD=DDD+1   !.. update day of year

         READ(1,'(13I3,3F5.1)',END=90,ERR=91) YEAR(K), MONTH(K),
     >   DAY(K),(AP(I),I=1,9),QB,F107ADJ(K),F107AVE,F107YR
      ENDDO
     
 90   CONTINUE
 
      !.. Now write out the results
      DO J=1,K-1
         !.. Some of the IRI average values are wrong. Recalculate the 81-day 
         !.. average F10.7 to replace IRI values.
         IF(J.GT.41.AND.J.LT.K-41) THEN
           CALL CALF107A(IDIM,J,F107OBS,F107ADJ,F107OBSB,F107ADJA)
           F107OBSA(J)=F107OBSB
         ENDIF
         
         !.. Output the new FLIP F10Kp values for the observed F10.7
         WRITE(3,'(I4,I3,I3,8I3,I4,I4)') YYYY(J),MONTH(J),DAY(J),
     >       (KP(J,I),I=1,8),NINT(F107OBS(J)),NINT(F107OBSA(J))
      ENDDO
      WRITE(3,92) 
 92   FORMAT(' ???? Check Last 13 lines of this F10KP-68.DDD for'
     >     ,1X,'bad Kp (e.g. all Kp = 0) ????')
      STOP
      CLOSE(3)

 91   CONTINUE
      J=K-1
      WRITE(3,*) ' '   !.. Inserts a blank line
      WRITE(3,*) '  **** ERROR Reading Input file near following date'
      WRITE(3,*) ' Year Month Day'
           WRITE(3,'(13I5)') YEAR(J),MONTH(J),DAY(J)
      CLOSE(3)
      STOP
      END
C::::::::::::::::::::: CALF107A :::::::::::::::::::::::::::::
C...  This function calculates the 81 day average F10.7 index
C...  Written by P. Richards, August 2008
      SUBROUTINE CALF107A(IDIM,K,F107OBS,F107ADJ,F107OBSA,F107ADJA)
      IMPLICIT NONE
      INTEGER I               !.. loop control variable
      INTEGER ISUM            !.. number of values in sum
      INTEGER K               !.. Array index for last value read in
      INTEGER IDIM            !.. Array dimension
      REAL F107OBS(IDIM)      !.. F10.7 index observed
      REAL F107ADJ(IDIM)      !.. F10.7 index observed
      REAL SUMOBS,SUMADJ      !.. F10.7 sums
      REAL F107OBSA,F107ADJA  !.. 81 day average F10.7s

      SUMOBS=0.0         !.. initialize sum    
      SUMADJ=0.0         !.. initialize sum    
      ISUM=0             !.. number of values in the sum

      !.. Sum the 81 indices
      DO I=K-40,K+40
        SUMOBS=SUMOBS+F107OBS(I)
        SUMADJ=SUMADJ+F107ADJ(I)
        ISUM=ISUM+1
      ENDDO

      !.. Calculate average F10.7
      F107OBSA=SUMOBS/REAL(ISUM)
      F107ADJA=SUMADJ/REAL(ISUM)

      RETURN
      END

C************************** SUBROUTINE CONV_DATE ***********************
C---- Scot A Braze,  8/7/95
C---- This subroutine converts the date from YYYY MM DD to YYYY DDD. 
C
      INTEGER FUNCTION CONV_DATE(YYYY,MM,DD)
      IMPLICIT NONE
      INTEGER MM,DD,NCONV(12),LYCONV(12),YYYY
      LOGICAL DEBUG
      DATA DEBUG /.FALSE./

      DATA NCONV/0,31,59,90,120,151,181,212,243,273,304,334/
      DATA LYCONV/0,31,60,91,121,152,182,213,244,274,305,335/

      IF(MOD(YYYY,4).EQ.0) THEN
          CONV_DATE=LYCONV(MM)+DD
      ELSE
          CONV_DATE=NCONV(MM)+DD
      ENDIF
      IF(DEBUG) WRITE(6,*) YYYY
      END

C:::::::::::::::::::::::: TEST_DATA ::::::::::::::::::::::::::::::::::::::::
C..  This subroutine tests to make sure that the year, month, and day, are 
C..  continuous
C..  Written by P. Richards, August 2008
      SUBROUTINE TEST_DATA(IDIM,K,YYYY,MONTH,DAY)
      IMPLICIT NONE
      INTEGER K           !.. Array index for last value read in
      INTEGER IDIM        !.. Array dimension
      !.. Year, Month, and day array
      INTEGER YYYY(IDIM), MONTH(IDIM), DAY(IDIM)

      !.. Test to see if there is a problem with the date
      !.. If year is the same, test Month then Day
      IF(YYYY(K)-YYYY(K-1).EQ.0) THEN 
        !.. Same month
        IF(MONTH(K)-MONTH(K-1).EQ.0) THEN
           !.. If days not contiguous
           IF(DAY(K)-DAY(K-1).NE.1) GOTO 20
        !.. Month has changed, test day=1
        ELSEIF(MONTH(K)-MONTH(K-1).EQ.1) THEN
           IF(DAY(K).NE.1) GOTO 20
        !.. Month difference not o or 1 - Something wrong
        ELSE
           GOTO 20
        ENDIF

      !.. Year has changed, test Month and Day
      ELSEIF(YYYY(K)-YYYY(K-1).EQ.1) THEN
        !.. Month should change from 12 to 1
        IF(MONTH(K-1).EQ.12.AND.MONTH(K).EQ.1) THEN
           !.. Day should change from 31 to 1
           IF(DAY(K-1)-DAY(K).NE.30) GOTO 20
        ELSE
           GOTO 20 !.. Month has changed by wrong amount
        ENDIF

      !.. Year has changed by too much
      ELSE
           GOTO 20 
      ENDIF

      RETURN    !.. If no errors return
20    CONTINUE
      REWIND(3)
      WRITE(3,192) YYYY(K),MONTH(K),DAY(K)
      CLOSE(3)
      REWIND(4)
      WRITE(4,192) YYYY(K),MONTH(K),DAY(K)
      CLOSE(4)
      WRITE(6,192) YYYY(K),MONTH(K),DAY(K)
      STOP

 192    FORMAT('  ** WARNING: Problem in Input file at YEAR=', I5,
     >       '  Month=',I3,'  Day=',I3,/
     >       '  Make sure that year, month, and day are continuous')     
      END
C:::::::::::::::::::::::: TEST_F107 ::::::::::::::::::::::::::::::::::::::::
C..  This subroutine tests to make sure that the F10.7 data are reasonable
C..  Written by P. Richards, August 2008
      SUBROUTINE TEST_F107(IDIM,K,YYYY,MONTH,DAY,F107ADJ)
      IMPLICIT NONE
      INTEGER K           !.. Array index for last value read in
      INTEGER IDIM        !.. Array dimension
      REAL F107ADJ(IDIM)  !.. Array of F10.7 values
      !.. Year, Month, and day array
      INTEGER YYYY(IDIM), MONTH(IDIM), DAY(IDIM)

      IF(F107ADJ(K).LT.50.OR.F107ADJ(K).GT.500) THEN
         REWIND(3)
         WRITE(3,192) YYYY(K),MONTH(K),DAY(K)
         CLOSE(3)
         REWIND(4)
         WRITE(4,192) YYYY(K),MONTH(K),DAY(K)
         CLOSE(4)
         WRITE(6,192) YYYY(K),MONTH(K),DAY(K)
         STOP
      ENDIF 
       
      RETURN
 192      FORMAT(//'  ** WARNING: Bad F10.7 index in input data file'
     >     /,4X,'at YEAR=',I5,'  Month=',I3,'  Day=',I3,/
     >     '    Edit file and replace bad value with the correct value.'
     >     /,4X,'OR, The previous day index may be reasonable'//)
       END
